setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM")
library(ggplot2)

result=read.table("20210318.Ref.Alt比例/117012.Ref.Alt.reads.txt",head=T,sep="\t")
affects=c("X1T","M7","M5","M1","M47","M49","M28","M27","M30","M29","M26","M25","M35","M36")
unaffects=c("X2B","M8","M6","M2","M48","M50","M18","M17","M20","M19","M22","M21","M40","M39")

sel1=c(paste0(affects,"_reads1"),paste0(affects,"_reads2"))
sel2=c(paste0(unaffects,"_reads1"),paste0(unaffects,"_reads2"))

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
library(openxlsx)
filea=read.csv("20201112做汇总表/all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210316LIBD.eQTL处理/53K.add.GWAS.eQTL.DEG.motif.for.analysis.txt",header=T,sep="\t")
file=file[!duplicated(file$id),]	#对TF去重
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

###可视化方面，考虑到数据中alt/ref计算方式将产生很多极值，箱型图整体会被极大值压缩，要使得图好看的话建议使用VAF画箱型图
##1.如果使用alt/ref，需去除部分极值，代码如下：
tmp=result[result$unitID %in% as.character(file3$unitID),]

data.affect=tmp[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$ratio.alt=data.affect$reads2/(data.affect$reads1)
data.affect=data.affect[is.finite(data.affect$ratio.alt),]
data.affect=data.frame(ratio.alt=data.affect$ratio.alt,group=rep("affect",dim(data.affect)[1]))


data.unaffect=tmp[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$ratio.alt=data.unaffect$reads2/(data.unaffect$reads1)
data.unaffect=data.unaffect[is.finite(data.unaffect$ratio.alt),]	###是否为有限值
data.unaffect=data.frame(ratio.alt=data.unaffect$ratio.alt,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data=data[data$ratio.alt<5,]	###去除部分极值

ggplot(data,aes(x = group, y = ratio.alt, fill = group)) +
  geom_boxplot(alpha=0.7) +
  scale_y_continuous(name = "ratio.alt")+
    scale_x_discrete(name = "") +
  ggtitle("200 ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 12),
        axis.text.x=element_text(size = 11))+scale_fill_lancet()
		
##################################################################################以下使用VAF画箱型图
##117012 ASH

data.affect=result[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$VAF=data.affect$reads2/(data.affect$reads1+data.affect$reads2)
data.affect=data.affect[is.finite(data.affect$VAF),]
data.affect=data.frame(VAF=data.affect$VAF,group=rep("affect",dim(data.affect)[1]))


data.unaffect=result[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$VAF=data.unaffect$reads2/(data.unaffect$reads1+data.unaffect$reads2)
data.unaffect=data.unaffect[is.finite(data.unaffect$VAF),]
data.unaffect=data.frame(VAF=data.unaffect$VAF,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data$group=factor(data$group,levels = c("unaffect","affect"))
p1=ggplot(data,aes(x = group, y = VAF, fill = group)) +
  stat_boxplot(geom = "errorbar",width=0.15)+geom_boxplot(aes(fill=group),outlier.colour="NA") +#加误差棒，加箱型图且设置离群值的颜色为空
  stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="white")+	#钻石表示均值
  scale_fill_manual(values = c("#00468B","#ED0000"))+
  scale_y_continuous(name = "VAF")+
  scale_x_discrete(name = "") +
  ggtitle("117012 ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 15),
        axis.text.x=element_text(size = 15))+guides(fill=F) 
##61725 ASH
tmp=result[result$unitID %in% as.character(filea1$unitID),]

data.affect=tmp[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$VAF=data.affect$reads2/(data.affect$reads1+data.affect$reads2)
data.affect=data.affect[is.finite(data.affect$VAF),]
data.affect=data.frame(VAF=data.affect$VAF,group=rep("affect",dim(data.affect)[1]))


data.unaffect=tmp[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$VAF=data.unaffect$reads2/(data.unaffect$reads1+data.unaffect$reads2)
data.unaffect=data.unaffect[is.finite(data.unaffect$VAF),]
data.unaffect=data.frame(VAF=data.unaffect$VAF,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data$group=factor(data$group,levels = c("unaffect","affect"))
p2 <- ggplot(data,aes(x = group, y = VAF, fill = group)) +
  stat_boxplot(geom = "errorbar",width=0.15)+geom_boxplot(aes(fill=group),outlier.colour="NA") +
  stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="white")+
  scale_fill_manual(values = c("#00468B","#ED0000"))+
  scale_y_continuous(name = "VAF")+
  scale_x_discrete(name = "") +
  ggtitle("61725 ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 15),
        axis.text.x=element_text(size = 15))+guides(fill=F) 
##53425 ASH
tmp=result[result$unitID %in% as.character(file$unitID),]

data.affect=tmp[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$VAF=data.affect$reads2/(data.affect$reads1+data.affect$reads2)
data.affect=data.affect[is.finite(data.affect$VAF),]
data.affect=data.frame(VAF=data.affect$VAF,group=rep("affect",dim(data.affect)[1]))


data.unaffect=tmp[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$VAF=data.unaffect$reads2/(data.unaffect$reads1+data.unaffect$reads2)
data.unaffect=data.unaffect[is.finite(data.unaffect$VAF),]
data.unaffect=data.frame(VAF=data.unaffect$VAF,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data$group=factor(data$group,levels = c("unaffect","affect"))
p3 <- ggplot(data,aes(x = group, y = VAF, fill = group)) +
  stat_boxplot(geom = "errorbar",width=0.15)+geom_boxplot(aes(fill=group),outlier.colour="NA") +
  stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="white")+
  scale_fill_manual(values = c("#00468B","#ED0000"))+
  scale_y_continuous(name = "VAF")+
  scale_x_discrete(name = "") +
  ggtitle("53425 ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 15),
        axis.text.x=element_text(size = 15))+guides(fill=F)
### 8544 ASH
tmp=result[result$unitID %in% as.character(file1$unitID),]

data.affect=tmp[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$VAF=data.affect$reads2/(data.affect$reads1+data.affect$reads2)
data.affect=data.affect[is.finite(data.affect$VAF),]
data.affect=data.frame(VAF=data.affect$VAF,group=rep("affect",dim(data.affect)[1]))


data.unaffect=tmp[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$VAF=data.unaffect$reads2/(data.unaffect$reads1+data.unaffect$reads2)
data.unaffect=data.unaffect[is.finite(data.unaffect$VAF),]
data.unaffect=data.frame(VAF=data.unaffect$VAF,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data$group=factor(data$group,levels = c("unaffect","affect"))
p4 <- ggplot(data,aes(x = group, y = VAF, fill = group)) +
  stat_boxplot(geom = "errorbar",width=0.15)+geom_boxplot(aes(fill=group),outlier.colour="NA") +
  stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="white")+
  scale_fill_manual(values = c("#00468B","#ED0000"))+
  scale_y_continuous(name = "VAF")+
  scale_x_discrete(name = "") +
  ggtitle("8544 ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 15),
        axis.text.x=element_text(size = 15))+guides(fill=F)
##807 psy-ASH
tmp=result[result$unitID %in% as.character(file2$unitID),]

data.affect=tmp[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$VAF=data.affect$reads2/(data.affect$reads1+data.affect$reads2)
data.affect=data.affect[is.finite(data.affect$VAF),]
data.affect=data.frame(VAF=data.affect$VAF,group=rep("affect",dim(data.affect)[1]))


data.unaffect=tmp[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$VAF=data.unaffect$reads2/(data.unaffect$reads1+data.unaffect$reads2)
data.unaffect=data.unaffect[is.finite(data.unaffect$VAF),]
data.unaffect=data.frame(VAF=data.unaffect$VAF,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data$group=factor(data$group,levels = c("unaffect","affect"))
p5 <- ggplot(data,aes(x = group, y = VAF, fill = group)) +
  stat_boxplot(geom = "errorbar",width=0.15)+geom_boxplot(aes(fill=group),outlier.colour="NA") +
  stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="white")+
  scale_fill_manual(values = c("#00468B","#ED0000"))+
  scale_y_continuous(name = "VAF")+
  scale_x_discrete(name = "") +
  ggtitle("807 psy-ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 15),
        axis.text.x=element_text(size = 15))+guides(fill=F)
##200 psy-ASH
tmp=result[result$unitID %in% as.character(file3$unitID),]

data.affect=tmp[,sel1]
data.affect$reads1=rowSums(data.affect[,paste0(affects,"_reads1")],na.rm = T)
data.affect$reads2=rowSums(data.affect[,paste0(affects,"_reads2")],na.rm = T)
data.affect$VAF=data.affect$reads2/(data.affect$reads1+data.affect$reads2)
data.affect=data.affect[is.finite(data.affect$VAF),]
data.affect=data.frame(VAF=data.affect$VAF,group=rep("affect",dim(data.affect)[1]))


data.unaffect=tmp[,sel2]
data.unaffect$reads1=rowSums(data.unaffect[,paste0(unaffects,"_reads1")],na.rm = T)
data.unaffect$reads2=rowSums(data.unaffect[,paste0(unaffects,"_reads2")],na.rm = T)
data.unaffect$VAF=data.unaffect$reads2/(data.unaffect$reads1+data.unaffect$reads2)
data.unaffect=data.unaffect[is.finite(data.unaffect$VAF),]
data.unaffect=data.frame(VAF=data.unaffect$VAF,group=rep("unaffect",dim(data.unaffect)[1]))


data=rbind(data.affect,data.unaffect)
data$group=factor(data$group,levels = c("unaffect","affect"))
p6 <- ggplot(data,aes(x = group, y = VAF, fill = group)) +
  stat_boxplot(geom = "errorbar",width=0.15)+geom_boxplot(aes(fill=group),outlier.colour="NA") +
  stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="white")+
  scale_fill_manual(values = c("#00468B","#ED0000"))+
  scale_y_continuous(name = "VAF")+
  scale_x_discrete(name = "") +
  ggtitle("200 psy-ASH") +
  theme_bw() +
  theme(plot.title = element_text(size = 15),
        text = element_text(size = 15),
        axis.text.x=element_text(size = 15))+guides(fill=F)
##
layout <- matrix(c(1,2,3,4,7,5,6,7), nrow = 2,byrow = T)
multiplot(plotlist=list(p1,p2,p3,p4,p5,p6),layout=layout)